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ROMBERG INTEGRATION 


Matthew J. O’Malley 
SUMMARY 

This paper presents a theoretical development of a pro- 
cedure to numerically integrate the definite integral 

fb 

I f(t) dt . Theorems and the majority of proofs are given, 

justifying the procedure, and remarks are made conceding the 
types of functions for which the procedure appears well suited. 

INTRODUCTION* 


This report presents a theoretical method to numerically 

fb 

integrate the definite integral I f(t) dt . A special 

•a 

case of the method, the Romberg Integration scheme, is also 
presented. Theorems and the majority of proofs are given 
justifying the procedure, and remarks are made concerning 
types of functions for which the procedure appears well suited. 

t 

Emphasis has been placed on the mathematical justification of 
the procedure in order to provide a deeper understanding of 

_ s 

the method, and, hopefully, to lead to further research of 
the procedure and its modifications. 
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the chief references from which the bulk of the material 
of this paper was obtained are Moler 1 and Bauer, Rutishauser, 
and Stiefel. 2 


The author wishes to thank Mr. Robert Meyers for the 
many helpful suggestions made during the preparation of this 
paper. 

SYMBOLS 


E 2n* n * 1,2 » 


E m» m * 0,1,** 

f n , n *= 0,1,** 

f Cn) , n * 1,2, 


“k 

M(h) 

N., k = 0,1, 


pj k) (x) 


p(k) 


1 


Bernoulli numbers 
real constants 
f (a + n h) 

n tn derivative of f 

integration step-size 

(b-a)/N k 

midpoint rule sum 

increasing sequence of positive 
integers 

(k-j) tk degree polynomial 
interpolating the points ^ 
(Xj, /*)» i - j > * * * »k 

value at x* of Pj k ^ (x) 
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value at 0 of the polynomial 
which interpolates the function 
x s ; (s > 1 , at the points 


' — 

Xj , i - 

T(0) 

r b 

J f(t) dt 

T(h) 

trapezoidal sum 

T 00 
• J 

value at 0 of (x) 

x A , i * 0,1, • • • 

distinct real or complex numbers 

y-, i >= 0,1,*** 

real or complex numbers not 


•N 

necessarily distinct 

a, p 

constants between 0 and 1 

S' 


t 





NEVILLE'S ALGORITHM 

Suppose^ that we are given (m+l) distinct points, 
x 0* x i » ***» x m » real or complex, and (m+l) corresponding 
values (not necessarily distinct), y 0 , y x , •••, y m . 
Neville's algorithm is a method of calculating the.- unique 
polynomial P(x) of degree m which takes on the values 
y k at the points x k . 

To describe the method, let pj k ^ (x) denote the inter- 
polating polynomial of degree (k-j) which satisfies 
Pj K) (x.) s , (i 15 j,**’,k) . 

(1.1) THEOREM. For j = k , let p£ k) (x) *= y k and 
for j < k , let 

Pj k) (x) = [(x k - X) pj k * 1 J (x) 

- (x^ - x) pj k J(x)j^(x k - Xj) (1.2) 

(kl 

P j (x) is then the unique polynomial of degree (k-j), 
which interpolates the points (x . , y^) , •••, (x k , y k ) . 

Proof. Induction is used on n « k - j . If n=0^ 
the theorem is obvious. Suppose the theorem has been proved 
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for n » m . Then 

- [C^ m ,i - x) 

- • *) p|^ nn) w]/( x j + mu - *J> 

• • l 

and by the induction hypothesis, pO +n O( x ) and (x) 

are polynomials of degree m , and hence Pp +m+1 )(x) is 
a polynomial of degree m+1 . Further, by direct verifi- 
cation, pjj +m+1 )( Xi ) = y^ for i * j, j+m+1 . By assump- 
tion, for j < i < j+m+1 , pO +m ) ( x .) = Pjli" 1 * 1 ^ (**) * y ± ; 
hence pj J+m+1 ^ (Xj) * y ^ for j < i < j+m+1 . 

Finally, the uniqueness follows from the fact that if 
P(x) and Q(x) are two polynomials of degree n interpo- 
lating the same C^x) distinct points, then P(x) - Q(x) is 
a polynomial of degree n having at least (n+1) roots and 
thus is the zero polynomial. Therefore P(x) « Q(x) . 

Neville’s algorithm then can be used to evaluate the 
interpolating polynomial P (x) at any desired point x* . 
An important advantage of Neville’s algorithm over the, per- 
haps, more familiar Lagrangian representation is that the 
number of points to be interpolated may be increased without 
redoing previous computation. For example, if we wish to 
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calculate pj 4 ^(x) , then it is only necessary to calculate 
pj 3 ) (x) and p[ 4 ) (x) and apply (1.2) to compute p£ 4 ^(x) . 
Thus, the degree interpolating polynomial p£ 4 ^(x) is 
determined'by using linear .nterpolation on ^x Q , P^ 3 ^(x)j 
and ^x 4 , p{ 4) (x)j . That this is true in general is clear. 
It is also important to realize that no assumption has been 
made regarding the distribution of the po’nts x^ . They 
need neither be equally spaced, nor in increasing or 
decreasing order. 


If we use the notation pj k ? for the computed va ues 

P? k) (x*) , then Neville's algorithm can be arranged in the 
3 

following table: 




p(°) 

r 0 

pCO 

r i 

pCD 

r o 


• 

p(2) 

r 2 

p(2) 

M 

p(2) 

r 0 


pO) 

*3 

p(3) 

p(3) 

M 

p( 3 ) 


• 

• 

• 

• 

• • 

• • 


• 

• 

• • 

p(k) 

F k 

■>) 

p (k) 

F k-2 

P (k) p (k) 
F k-3 F o 




( 1 . 3 ) 


■s 



is obtained directly 
immediately to the 


Using this arrangement each entry Pp) 

from the two entries pP" 1 ) am. Pp2 

1 1 + 1 

*'leftl!L-and fi left“ above M it. 


In the usual case, T(x) is a function such that 
T(x^) = y^ for C < i < n . To approximate T(x*) by using 
polynomial interpolation, we iteratively compute the values 
Pp) until two or more successive values, P^) t pp + l) , 

*•* t agree to within some preassigned degree of accuracy. 

Finally, we note that if x* . « 0 , 




p(k-l) 

p i 



( 1 . 4 ) 


We shall be interested in‘ this form when using the Romberg 
scheme . ' ■ 


* 


s 



ROMBERG INTEGRATION 


One of the most well-known methods of approximating 

fb 

the integral J f(t) dt is the trapezoid rule. If N is 


"a 


a positive integer and 


h * fb - a)/N 


x n ■ a + nh 


f n “ £ ^ x n^ » n * 0,1,2, ••• ,N 


then the N-interval trapezoid rule determined by the 
N-subintervals , (x^, x^ +1 ) , i * 0,1,*»*,N-1 , is given by 


T(h) * f„ ♦ fj + 


+ '»■! * J f N ) 


( 2 . 1 ) 


where h is the mesh size. 

r b 

It follows easily, then, that if j f(t) dt exists, 

b ** 

Tfh) f f(t) dt as h -*■ 0 . Further, if f" (t) is con- 
•'a 

tinuous on [a,b] , and hence bounded on [a,b] , then 


T(h) = f b f(t) at + l(b - a) h 2 /12] • f"(S) 


where C e (a,b) . 



The Romberg method consists of an application of 
Neville's algorithm to the function T(h) . In general, let 
N p , N k , ••• be an increasing sequence of positive 

integers and 

h k = (b - a)/N k 

x k * h k (2.2) 

• ' Y k = T(h k ) , k = 0,1,..- 


fkl 

Applying Neville's algorithm with x* * 0 and TV 1 denoting 
the value P^) , we hcve from (1.4), 

For j = k , T^ T(h k ) and for j < k , 

s 

T j k) ■ [ h k T f‘ i; - h j T jn]/( h k - h j) 


(2.3) 


Thus, TV^ is the value at 0* of the (k-j)** 1 degree poly- 
nomial which interpolates the (k-j+1) points, fhj, T(h^)j , 
•••,* |h k , T(h k )j. . Therefore, from (1.3) we have the 


following table: 



\ 




T (l) 

T (2) 

*1 

T (3) 
1 2 



T (2) 

*0 

t (3) 



t (3) 

l C 

•' • 

• • - 
• • 

T (k) T (k) 
*k-3 *0 



( 2 . 4 ) 


• • • * • . ♦ •_ 

• • • • • 

Each entry in the table is, an approximation to T(0) ; 

f b 

that is, I £(t) dt . The first column (which we shall 
J a 

Call the zero tb column) contains the values of the succes- 
sive trapezoid rules. We should note here that the process 
just described is often called Romberg Integration. However, 

we shall restrict ourselves to. this terminology for only the 

v 

special case - 2 . The general procedure just described 

will be referred to as a modification or generalization of 

the Romberg scheme. As we shall see in Section III, the 

basis for this procedure will be the existence of an asymp- 

2 * 

totic expansion of T(h) in powers of h . For this 
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reason, the above process is often referred to as extrapo- 
lation to the limit or Richardson’s deferred approach to 
the~limit. 3 

As noted previously, the Romberg scheme uses the values 

Ir 

N^-2 . This choice has several advantages 

for computational purposes. For h * (b - a)/N , let 
M(h} be the N- interval midpoint rule; that is 

MCh) » h f(a ♦ (" - \) h) 

Since 

* » 

T(h) = ffa * n h) ♦ (1/2) (f(a) * f(b))J 


it follows that T(h/2) = (T(h) + M(h))/2 . For the Romberg 
scheme, since h^ +1 = h^/2 , it follows that 


T(W (T(h k ) *M(h k ))/2 


(2.5) 


This relation was used by Romberg to construct the zero 
columr of his array. Furthermore, using the fact that 
hj/h^ = 2^*3 , (2.3) can be written in the form: 


th 


T (k) 


T j+1 + ( T j+i ' T j kl) )/( 4k ’ J ’ 1) (2-6) 
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In jthis form, th e alg orithm was first described by Romberg 
^ — and is called Romberg Integration. 4 

1c 

However, while the choice * 2 has many "nice" 
computational features, it also has the disadvantage that the 
number of grid points at which f must be evaluated doubles 
with each iteration. If f is a complicated function, this 
could be a significant disadvantage. 

The natural suggestion would be to choose a sequence 

( v)~ 

N v less rapidly increasing than the sequence {2 } 

K { Jk*=0 

However, before we consider such choices, we should concern 

ourselves with the question of when does the Romberg scheme 

r b 

converge to I f(t) dt . We shall consider this problem 
’'a 

and others in the next section. 
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CONVERGENCE CRITERIA 


Our first aim in this section will be to justify the 
Romberg scheme. We will show that every diagonal and every 
column of (2.4) with N^ = 2^ (k = 0,1,**») converges to 
T(0) . Proofs are based on those in ref. 2. 

It is clear that every entry of (2.4) is a linear com- 
bination of elements of the zero^ 1 column; that is. 


T (k+m) _ v^m T (k+i) 

T k “ i=0 c m,m-i l k+i 


(3.1) 


where the coefficients c m m _-. are independent of k . 
From (2.6), it follows that for m j* 0 , 


T (k+m) 

*k 



T (k*m) 

x k+l 


T Ck.m- 1 ))/ (4 m . 1} 


(5.2) 


From this observation and from (3.1), it follows easily by 
induction on m that for m f 0 , the coefficients c m m _ - 
obey the recursion formula 


; . « [ 4 “ c , . - c . . -l/(4 m - 1) (3.3) 

m,m-i l m-l,m-i m-l,m-l-ij/ v 


for i = 0,1, . assuming c H . l m - ■= 0 . 
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If t m (x) are the polynomials defined by 
. ^”-0 <=»,* * k • --0,1,2,... (3.4) 

then from (3.3) we have that, for m f 0 , 

* • - t 

* b M ' Ek.op 4 * c m-l,k " Si-l.k-l)/^” * l) j xk 

- [(4 m - x)/(4 m - 1)] Vl k x k 

- (4 m - X) t B .j(x)/(4 m - 1) 

✓ 

s 

From this result it follows that, for m ^ 0 , 

t B (x) = - xJ ]/[ n i-l C4 1 - 1)] (5.5) 

which, together with the fact that c Qj0 c 1 * a H°ws us to 
conclude that 


V 1 ) 


E m 

i=0 


m, i 


V a c 

m,m-i 

V. ' 


1 


(3^6) 


for each m . 


is 


Further, we have the following: 

(3.7) LEMMA. For each m , Ei=o| c m i| < 2 
Proofs' Since 

= [n? =1 (4 i - u]/[n? =1 C4 1 - 

£ [n" =1 (4 i ♦ 0/(4* - i)J < 2 

it suffices to show that Y*? n |c . * t (-1) for each m 

" i=u I i.ii m v J 

. ' • 

To prove this, it is sufficient to show that (-1) 1 . > 0 

for i = 0,1,2, . Suppose that m = 1 . Since 

(-1)° c n - c n * t (0) > 0 for all m , we have that 

(-1)° c lj0 > 0 . Further, using (3.1) and (3.2), it can be 

observed that c j ^ = -1/3 ; hence (-1) 1 c^j * 1/3 > 0 . 

Suppose the statement has been proved for m = n - 1 . We 

consider (-1) 1 c . , where i e {0,1, • • • ,n-l} , From 

n, l 


(3.3), 

<-»* c n,i ' c n-l,i 



c n-l,i-l 



and by the induction hypothesis, (-1) 1 c n _i \ > ® > an< * ** 
(-1)** c n _i > 0 i e {0, !,•••, n-1) . Therefore, 
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[«"(-« 1 c n-l,i * C - 1 ) 1 ' 1 Vl,i-l| (4 ’ - 1 ) > 0 . Finally, 
r S-l.« c „-l,n-l j/t 4 ” • « 

- [o ♦ C-l )"' 1 C n .i >n .i]/c 4 n - 1) > 0 


c-l) n C. 


,n-l 


since (-1) c n -l n-1 > P by the induction hypothesis. 

•Before the convergence of the Romberg method can be 
shown, we need to observe that for each k , 


lim c. 
m-»“> 


in,m-k 


(3.8) 


. To show this, observe that from (3.5) 


t m C4x) - 4 m t m (x)(l - *)/(4" - x) 




and 


Hence 


or 


Vl (x) * t m O°(4 m - D/C 4 ” - x) 


V 4x > - V x) * ' x Vl (x) 


Ek-o ' « xk ' -SwW^ 1 
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By equating coefficients, we then have 

c m,k " " c m-l,k-l " 1) 

Repeated application of this last relation gives 

■ Vm-k - c k>0 /n?:i (4 1 - 1) (3.9) 

From this last relation, (3.8) is apparent. 

The convergence of the Romberg scheme may now be 
proved. 


(3.10) THEOREM. The convergence of the zero** 1 column 
of the array (2.4) implies the convergence of all further 
columns of (2.4) to the same limit* 

/* 

Proof. Suppose that lim t£*^ = T(0) = f f(t) dt . . 

k+® X J a 


Consider the elements , k = 0,1,* •• 

column where •pCk+m) V'm _ 'p(k'^i) 

column, wnere i k 2v i=0 c m,m-i T k+i 

We show that lim = T(0) . Let e > 0 

k-*°° K 


, of the m th 
for each k . 
be assigned. 


Choose N , a positive integer, such that 
| T 00 _ T(0) | < e/ [ (m ’+ 1) U] for k >N , where 

U ' max | c m,m-i| • 1 = • By (3-6) 


s 


"nijm-i 


i 
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hence 

T (k+m) . T ( 0) 


Therefore, for k > N , 

| T (k+m) . T ( 0 )| < (m + u • e/[(m + 1) • U] * e 

and the theorem is proved. 

(3.11) THEOREM. The convergence of the zero 1 * 1 column 
of the array (2.4) implies the convergence of all diagonal 
sequences; that is, of all sequences T^ k+m ^ , k constant, 
m <° , to the same limit. 

Proof. Let k be a fixed nonnegative integer. Under 

the assumption that lim T^ = T(0) = f f(t) dt , we 

k-»-» K 

show that lim T^ k - T(0) . We note that a sequence- to.- 
sequence transformation is established by (3.1) between the 


'2". 0 T k*i i5 ' T C°) 


XU • s;.0 V-i rm 

£i"0 c m,m-i( T k.t i) - T (0)) 
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sequences (t^L - 

to the Silverman-Toeplitz theorem (Theorem 4.1, II, p. 64), 5 
a necessary and sufficient condition that x(^ +m ) -► j(0) , 
as m «> , is that 


(a) 

£i*0 | c m,m-i| 

< M 

for 

every positive integer 

(b) 

r 

lim c . = 0 

m-*» ®» m "h 

for 

each 

fixed k , and 

(c) 

V* n c 
4'i=0 m,m - 1 

A -»■ 
m 

1, 

as m -*■ » . 


These conditions are satisfied by Theorem (3.7) and statements 
(3.6) and (3.8). 


.. Then, according 

m=0 


Thus, under only the assumption that f is Riemann 
integrable on [a,b] , the convergence of the Romberg scheme 
has been shown. Therefore, in the usual case of interest, 

noted in Section I, the diagonal entries T^) converge to 

fb 

T(0) if only I f(t) dt exists. However, it has not yet 
•'a 

been shown which column or diagonal converges to T(0) 
"faster.’’ The remainder of this section will be concerned 
with this problem. 

v S 

S 

* / * 

/ 
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(3.12) Definition. A function f(x) has "order 

v 

as x approaches 0," and we write f(x) - 0(x ) , if con- 
stants M and x Q exist such that for all |x| < x 0 , 

|f(x)| £ M|x k | . 

(3.13) Definition. A function f(x) has an asymptotic 
series as x approaches 0, if constants a Q , , ••• exist 
such that for all m 

y(x) = a 0 + a x x + a 2 x 2 + • • • + a ffi x m + 0(x 111 ) (3.14) 

(3.15) THEOREM. If fC 2m+2 )( t ) is continuous for 
a <t < b , then there are numbers a 2 , •••, a m depending 

on f , but not on h , so that 

T(h) = J b f(t) dt + a t h 2 + ••• + a ffl h 2m + 0(h 2m+2 ) 

(3.16) 

In other words, T(h) as an expansion consisting of the 
first (m+1) terms of an asymptotic series. 

Proofs of this result will be found in Ralston 6 and 
9 7 

Henrici (Theorem 13.6, p. 257). We can now prove a theorem 

which is of a more general nature than the Romberg scheme and 

which will describe the "speed of convergence" of the columns 

fkl 

of the array (1.3), We will consider the quantities PV J 
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as j and k both tend to infinity, with n = k - j held 
fixe'’ These form the n** 1 column of the array (1.3), con- 
tinuing to designate the first column as the tero 1 ^ column. 
The theorenTand the proof given here are in ref. 1. It is 
essentially the same theorem as that given in ref. 7 
(Theorem 12.4, p. 240), except that in the latter theorem, 
it is assumed that the ratio (hj^ is a fixed constant 
between 0 and 1 for all k . 


(3.17) THEOREM. Assume that for m , a positive 
integer, the function y(x) has^an expansion consisting of 
the first m+1 terms of an asymptotic series; that is, 

y(x) = a Q + aj x + ••• + a ffl x m + 0(x m+1 ) (3.18) 


as x approaches 0. Assume further that the evaluation 
pcints x Q , x x , ••• satisfy 

0 < o £ x k/ x k-l < P < 1 (3.19) 

for all k . Then, for n < m , the values p£^ defined 
by Pj^ * ^k-n^ 0 ^ with j = k - n satisfy 


P 00 » 
k-n 


a. 


(-D n 


a n+l x k-n 



(3.20) 
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as k approaches 00 . Furthermore, if a n i 0 , then 

k([-k ■ -.]/[«.. • •.]) - • 


That is , the n** 1 column converges to a 0 faster than the 
(n-l) st column. 

Moler 1 proves the theorem by first proving the following 
lemma. 

(3.22) LEMMA. Assume x Q , x 1 , ••• satisfy (3.19) and 
for any positive integer s define QV ' by 

O 00 . x s 

Q k,s x k 

and for j < k , - ' / 

<; • Mr - >, AX.]h ■ 

Then 


x k-n,: 



. if 

s 

< 

n 

(3.23) 

x k-n "• V if 

s 

s 

n + i 

(3.24) 

as k •> » , if 

s 

> 

n + 2 

(3 v2$) 



is 

I 

\ . • 7 

Proof. By Theorem (1.1) Q^-n s * s t ^ ie va ^ ue at ® of 
the polynomial of degree n which interpolates the function 
x s —at— the points x^^, •••» x^ . 

Case 3. s £ n . s ( x ) " x S is a polynomial of 

degree n , which has n+1 zeros. Hence, it is identically 
zero and therefore Qj^ s * 0 . (In fact, s '^~n .) 

Case 2. s = n + 1 . cC x ) * x n+ * is. a polynomial « 

of degree n+1 . which has the n+1 roots, x^^, x^ . 

Thus, by the fundamental theorem of algebra, 

«*££,,« • * n+1 ■ - ** - w ••• ( X - x k ) 


Therefore, 


qCW 

x k-n,s 


- l(-l ) n+1 X 


• • • • Xm * 

k-n : x k 


(-l) n Xk_ n ... 


Case 3. s > n + 2 . We use (3.19) and induction on 
n . Let n = 0 . Then Q^ k g = x^ * O^x^ j . Assume the 
statement has been proved for n = m - 1 . Then 


pi 10 

x k-ra,s 



(k- 

k-m, 



s 
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By the induction hypothesis, we have 


- 

x k-m,s 

[ x k °( x k-i’ 

) - x k-m °( x k)]/ fx k ‘ x k- m > 


[('kK-m 1 

°( x k-l) - °( x k)]/( (x k/ x k-m ) * 

Since 0 < < 

1 £ x k/ x k-: 

l = °( x k) * Further » sil 

x k/ x k-l 5 p 

< 1 , we 

have that x j c / x i t _ ni ^ P™ an< * hence 

Q<» 

H k-m,s 

* [(1 + 

p")/Cl - ?")] • o(x*) = o(x») 


Therefore, the proof of Lemma (3.22) is complete. 


Proof of Theorem (3.17). We first establish the 
following proposition which is stronger than (3.20). Under 
the hypothesis of Theorem (3.17) 


k-n 


,0c) 


a + a , Qi J + • 
o n+1 y k-n,n+l 


+ a 


m 


x k-n,m 


+ 0 


(* 


m+l\ 

k I 


(3.26) 


Proof is again by induction on n . For n = 0 , the 
statement follows from (3.18) and the fact that Qj^ = 
Assume statement (3.26) is true for the case n-1 . Then.. 




.(k- 

k- 


k k-n 


P (k) 
*k-n+l 


h 


L k-n J 
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and by the induction hypothesis, we have 


p 00 

F k-n 


hi 




(k-1) , 

, v J + * • * 4 


a 0 (kl) + o/x m+1 \l 
a m "k-n,m + °\ x k-l/J 


t-nf 


x k-n| a o + a n ^k-n+l,n 


Ck) 


+ . 


♦ o^r 1 )])/(* 


k - W 




a. + a qP^ + a . , + * * • + a 

o n v k-n,n n+1 y k-n,n+l m y k-n, 

* [<*k/*k-n> °(*El) * °( x r 1 )l/(( x k/ x k-n ) - l) 


Since n * 0 by (3.23), and by the proof of (3.25), 

statement (3.26) is true. Finally, since 

i(n) 


and 


Qk-n,n + l = C-D“ x k-n 


«£Ui - °(*r) 


for i ** 2,***,m-n , (3.20) is true. 



Finally, if a n f 0 , then 
[ p k k i^ a o]/[ p k k Ll • a o] 

* [c- 1 ^ Vl x ’-n "• x k * °( x k +2 )]/ 

[(-l) n+1 a„ x k . n+1 ••• * k + °( x k* 1 )] ~ 

' [• 'VllV x k-n * °( , ‘)]/( } * ° tS k>) 

= Ofa^) » which implies (3.21). 

Therefore, if a n f 0 , the n** 1 column converges faster 

c + 

than the (n-1) column. 

This theorem is now applied to the general scheme 
developed in the Romberg Integration Section. Using the 
fact that the coefficients a n occurring in the expansion 
of T(h) in (3.16) have the form 

[(B 2n )/[(2n) l]j[f C2n ' 1) (b) - f< 2n -n (a) ] 
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where the B 2 n are the Bernoulli numbers, we have 

(3.27) COROLLARY. If f ^ 2m+2 ^ (t) exists ad is con- 
tinuous for a ^ t < b , and if 0 < a £ ^ P < ’ 

for all -lT, then for n < m , the entries in the n tb 
column of (2.4) satisfy 

T k-i - j[ bf(t)dt 

4 <-«" Ww— V 2 ♦ °( h k nt4 ) as k * * 

. (3.28) 

Furthermore, if 

f(2n- 1 ) ( . b ) ^ f (2n-l) (a) (3.29) 

th r b 

the entries in the n l column converge to I f(t) dt 

st ® 

faster than those in the (n-1) column. 

If we further assume that f is analytic on [a,b] , 

8 

then it follows from results of Gragg and Bulirsch and 
9 

Stoer that 

(3.30) THEOREM. If T(h) has an asymptotic expansion 

in powers of h (guaranteed if f is analytic on [a,b]), 

and if 0 < a £ /^k- 1 ^ P < 1 ' for all k , then for • 

each m > 0 , constants E„ exist such that 
* m 

■ l T « (k) - a »i S E B*l< h k •" Vm’ 2 < 3 ‘ 51 > 



Hence, if (3.28) and (3.31) are true, and » Vl * 0 • 
then the principal diagonal converges to a Q faster than 
the m^ ^-column of the array (2.4). Therefore, if a ffi f 0 
for all m > 1 , then the principal diagonal converges to 
a Q faster than any column of the array (2.4). 

In ref. 1 it is noted that Corollary (3.27) suggests * 
which functions may not work well with the scheme. For 
example, functions whose low- order derivatives do not 

* 

exist at the end points of integration, would make the 
expression (3.29) for the corresponding low-order coefficients 
meaningless. Also, periodic functions whose odd-order 
derivatives are equal at a and b we might expect not to 
be adaptable to the scheme. On the other hand, functions 
which are analytic or have high-order derivatives on [a,b] , 

t 

and are not periodic, we would expect to be, in many cases, 
particularly adaptable to this scheme. 

Bauer, Rutishauser, and Stiefel 2 remark that the argu- 
ments used to prove Theorems (3.10) and (3.11) can be 

h k-l £ 0 < 1 

for all k . Further, it is noted in ref. 2 that Jemma 
(3.7) guarantees the numerical stability of the Romberg N 
scheme, but for more slowly increasing sequences , the 

susceptibility to round-off error is increased. 


applied to the more general case if only h^j 


CONCLUDING 'REMARKS 


As noted previously at the end of the Romberg Integra- 
tion Section, it would seem practical to choose a sequence 

V 

less rapidly increasing than the sequence 2 . One 

choice that would certainly minimize the number of calcu- 
lations at each step would be *= k + 1 , k = 0,1,*** . 
However, this choice obviously does not satisfy the require 
ment £ p < 1 for all k and, as noted in ref. 2, 

can cause severe numerical instability if many columns of 
the table (2.4) are computed. 

Another choice, suggested by Bulirsch 10 is 

! 1 k = 0 

2 (k+l)/2 # k odd * . 

3 • 2 k/2_1 , k even 

That is, = 1, 2 ,3,4 ,6 ,8,12 ,16 , • • • . In this case 
P - 3/4 . 

Another choice suggested in ref. 2 is the sequence 
N^ = 1 ,2 ,3,6,9 ,18 ,27, 54, •• • (all powers of 3 and their s 

k-1 


doubles) . For this choice 


< p < 1 for all k , 
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but the modification is more susceptible to round-off errors 
than the Romberg scheme, since the sum ]T)?_q | c^ m _.| is 
higher than that for the Romberg method (up to 3.5). . 

The best choice of sequences appears to depend on the 
function, or at least class of functions, being integrated. 
In ^ref . 2 the function 1/x is integrated between 1 and 2 
using the 1,2, 3, 6, 9,* • • modification. It is noted that 
even if no improvement on the Romberg scheme can' be made, it 
does appear that computing time will be reduced to some 
degree in the long run. 

There is a very extensive derivation and discussion of 
the Romberg algorithm and its generalizations in ref. 2. 
Bulirsch 10 has extended the method’ of Romberg to any 
sequence of numbers h. with lim h. = 0 , by considering 

l"* 00 . 

a linear matrix transformation of the sequence T(h^) . 

Meir and Sharma 11 give an improvement of the result of 
Bulirsch. Stroud 1 ' has considered error estimates for the 
Romberg procedure, and compared them with error estimates of 
the Gaussian formulas. In a recent paper by Bulirsch and 

t 

1 3 

Stoer, it is proposed to use a Romberg- like scheme based 
on rational function extrapolation. The authors give 
applications of their method in refs. 13 and 14, and compare 



it with the Romberg scheme using polynomial interpolation in 
refs. 13 and 15. 

___The Romberg algorithm and its generalizations used for 

the numerical integration of definite integrals are based 

on the assumption that the trapezoidal approximation with 

step h has an asymptotic expansion in powers of h . 

It is proposed in refs. 2 and 9 to apply similar ideas to 

the solution of first-order ordinary initial-value problems 

using Euler's method as the basic discretization. The 

corresponding asymptotic expansion then also contains odd 

powers of h . Gragg 16 has established the existence of 

simple discretizations of both first and special second-order 

2 

systems which have asymptotic expansions in powers of h 

He then proposes to apply the modification of Bulirsch and 
< 

Stoer to obtain the solution for this type of ordinary 
initial-value problem. Refs. 2 and 1 give examples where 
the basic assumption of the existence of ?n asymptotic 
expansion of T(h) in powers of h is not valid. 



REFERENCES 


1. Moler,,Cleve B.: Extrapolation to the Limit. Numerical 

Analysis Seminar, University of Michigan Engineering 
Summer Conferences, June 19-30, 1967. 

2. Bauer, F. L. ; Rutishauser, H.; and Stiefel, E.: New 

Aspects in Numerical Quadrature. Proceedings of 
Symposia in Applied Mathematics, vol. 15, American 
Mathematical Society, 1963, pp. 199-218. 

3. Richardson, L. F. : The Deferred Approach to the Limit, 

I-Single Lattice. Philos. Trans. Roy. Soc. London. 

Ser. A, vol. 226, 1927, pp. 299-349. 

4. Romberg, W. : Vereinfachte Numerische Integration. 

Norske Vid. Selsk. Forn. (Trondheim), vol. 28, 1955, 
pp. 30-36. 

5. Cooke, Richard G. : Infinite Matrices and Sequence 

Spaces. Dover Publications, Inc., 1955. 

r 

6. Ralston, A.: A First Course in Numerical Analysis. 

McGraw-Hill, 1965. 

7. Henrici, P.: Elements of Numerical Analysis. John Wiley 

and Sons, Inc., 1964. 

8. Gragg, W. B.: Repeated Extrapolation to the Limit in 

the Numerical Solution of Ordinary Differential 
Equations. Doctoral dissertation. University of 
California, Los Angeles, 1964. . 

9. Bulirsch, R. ; and Stoer, J. : Fehlerabschatzungen ur.d 

Extrapolation mit Rationalen Funktionen bei Verfahren 
vom Richardson-Typus . Numer. Math., vol. 6, 1964, 
pp. 413-427. 



' ‘ • I : 

' # 33 

- - ' ! 

\ 

10. Bulirsch, R. : Bemerkungen zur Romberg- Integration. 

Numer. Math., voi. 6, 1964, pp. 6-16. 

11. Me ir , A. ; and Sharma, A.: On the Method of Romberg 

Quadrature. J. Siam Numer. Anal. Ser. B, vol. 2, 
no. 2, 1965, pp. 250-258. 

12. Stroud, A. H. : Error Estimates for Romberg Quadrature. 

J. Siam Numer. Anal. Ser. B, vol. 2, no. 3, 1965, 

pp. 480-488. ~ 

13 . Bulirsch, R. ; and Stoer, J.: Numerical Treatment of 

Ordinary Differential Equations by Extrapolation 
Methods.* Numer. Math., vol. 8, 1966, pp. 1-13. 

14. Bulirsch, R. ; and Stoer, J.: Asymptotic Upper and 

Lower Bounds for Results of Extrapolation Methods. 
Numerische Mathematik, vol, 8, 1966, pp. 93-104. 

15. Bulirsch, R. ; and Stoer, J.: Numerical Quadrature by 

Extrapolation. Numerische Mathematik, vol. 9, 1967, 
pp. 271-278. 

16. Gragg, W.: On Extrapolation Algorithms for Ordinary 

Initial-Value Pioblems. J. Siam Numer. Anal. Ser. B, 
vol. 2, no. 3, 1965, pp. 384-403. 


